library(tidyverse)
library(wesanderson)
library(MetBrewer)
library(sp)
library(sf)
library(classInt)
library(raster)
library(exactextractr)
library(fixest)
library(rgdal)
library(raster)
library(classInt)

# map of wealth, pm, scatter in corner, and 


############################################# functions

#add transparency to any color	
add.alpha <- function(col, alpha=1){
  if(missing(col))
    stop("Please provide a vector of colours.")
  apply(sapply(col, col2rgb)/255, 2, 
        function(x) 
          rgb(x[1], x[2], x[3], alpha=alpha))  
}	


#########################
get_grid_raster <- function(raster){
  output <- raster(raster)
  output[] <- 1:ncell(output)
  return(output)
}



identify_non_cells <- function(raster, poly){
  raster = raster(raster)
  raster[] <- 1:ncell(raster)
  r_ext <- extract(raster, poly)
  all_cells <- 1:ncell(raster)
  non_cells <- all_cells[all_cells %in% unlist(r_ext) == F]
  return(non_cells)
  
}
#########################

#############################################

    #load shapefiles

        #Rio, Katmandu, Lagos
        
        
        rio <- readOGR("data/city-boundaries/gadm40_BRA_2.shp")
        rio <- rio[rio$NAME_2=="Rio de Janeiro",]


        lagos <- readOGR("data/city-boundaries/gadm40_NGA_1.shp")
        lagos <- lagos[lagos$NAME_1=="Lagos",]


        delhi <- readOGR("data/city-boundaries/gadm41_IND_3.shp")
        delhi <- delhi[delhi$NAME_3=="Delhi",]

  
        
      #next section of code takes raw inputs and writes out data used for making figures.
      #input data are publicly available and can be downloaded for this section to run
        
        # ###########################get PM ######################
        # 
        # grid <- raster(raster("1_input/data/vandonk/annual_001/V5GL02.HybridPM25.Global.199801-199812.nc"))
        # grid[] <- 1:ncell(grid)
        # 
        # fls <- list.files("1_input/data/vandonk/annual_001/")[18:23]
        # rdat <- raster::brick(x = grid, nl = length(fls))
        # 
        # for (i in 1:length(fls)){
        #   
        #   rdat[[i]] <- (raster(paste("1_input/data/vandonk/annual_001/", fls[i], sep = "")))
        #   
        # }
        # 
        # 
        # pm_rio <- raster::brick(x = raster::crop(grid,extent(rio)), nl = length(fls))
        # pm_lag <- raster::brick(x = raster::crop(grid,extent(lagos)), nl = length(fls))
        # pm_del <-   raster::brick(x = raster::crop(grid,extent(delhi)), nl = length(fls))
        # 
        # for (i in 1:length(fls)){
        #   
        #   pm_rio[[i]] <- raster::crop(rdat[[i]], extent(rio))
        #   pm_lag[[i]] <- raster::crop(rdat[[i]], extent(lagos))
        #   
        #   pm_del[[i]] <- raster::crop(rdat[[i]], extent(delhi))
        # }
        # 
        # 
        # pm_rio <- raster::stackApply(x=pm_rio,indices = 1, fun= mean)
        # pm_lag <- raster::stackApply(x=pm_lag,indices = 1, fun= mean)
        # pm_del <- raster::stackApply(x=pm_del,indices = 1, fun= mean)
        # 
        # 
        # 
        # 
        # pm_rio_grid <- get_grid_raster(pm_rio)
        # pm_lag_grid <- get_grid_raster(pm_lag)
        # pm_del_grid <- get_grid_raster(pm_del)
        # 
        # noncell_rio <- identify_non_cells(pm_rio_grid, rio)
        # pm_rio[noncell_rio] <- NA
        # 
        # noncell_lag <- identify_non_cells(pm_lag_grid, lagos)
        # pm_lag[noncell_lag] <- NA
        # 
        # noncell_del <- identify_non_cells(pm_del_grid, delhi)
        # pm_del[noncell_del] <- NA
        # 
        # 
        # 
        # 
        # 
        # 
        # ###########################get wealth ######################
        # 
        # w_rio <- read_csv("1_input/data/rwi/BRA_relative_wealth_index.csv") %>% arrange(longitude, latitude)
        # bra_loc <- SpatialPoints(w_rio[,c("longitude","latitude")])
        # crs(bra_loc)<-crs(rio)
        # bra_loc$rwi <- w_rio$rwi
        # overlap <- over(bra_loc, rio)
        # w_rio <- w_rio[!is.na(overlap[,1]),]    
        # 
        # 
        # 
       # 
        # w_lag <- read_csv("1_input/data/rwi/NGA_relative_wealth_index.csv") %>% arrange(longitude, latitude)
        # lag_loc <- SpatialPoints(w_lag[,c("longitude","latitude")])
        # crs(lag_loc)<-crs(lagos)
        # lag_loc$rwi <- w_lag$rwi
        # overlap <- over(lag_loc, lagos)
        # w_lag <- w_lag[!is.na(overlap[,1]),]    
        # 
        # 
        # 
        # w_del <- read_csv("1_input/data/rwi/IND_PAK_relative_wealth_index.csv") %>% arrange(longitude, latitude)
        # del_loc <- SpatialPoints(w_del[,c("longitude","latitude")])
        # crs(del_loc)<-crs(delhi)
        # del_loc$rwi <- w_del$rwi
        # overlap <- over(del_loc, delhi)
        # w_del <- w_del[!is.na(overlap[,1]),]    
        # w_rio = w_rio
        # 
        # pm_rio = pm_rio
        # pm_lag = pm_lag
        # pm_del = pm_del
        # 
        # ###########################elevation ######################
        # 
        # el_rio1 <- raster("1_input/data/elevation/rio/ASTGTMV003_S23W044_dem.tif")
        # el_rio2 <- raster("1_input/data/elevation/rio/ASTGTMV003_S24W044_dem.tif")
        # el_rio <- raster::merge(el_rio1, el_rio2)
        # el_rio <- raster::crop(el_rio, extent(rio))
        # rm(el_rio1, el_rio2)
        # not_cell_rio <- identify_non_cells(el_rio, rio)
        # el_rio[not_cell_rio]<-NA
        # 
        # el_lag1 <- raster("1_input/data/elevation/lagos/ASTGTMV003_N06E002_dem.tif")
        # el_lag2 <- raster("1_input/data/elevation/lagos/ASTGTMV003_N06E003_dem.tif")
        # el_lag3 <- raster("1_input/data/elevation/lagos/ASTGTMV003_N06E004_dem.tif")
        # el_lag <- raster::merge(el_lag1, el_lag2)
        # el_lag <- raster::merge(el_lag, el_lag3)
        # el_lag <- raster::crop(el_lag, extent(lagos))
        # rm(el_lag1, el_lag2, el_lag3)
        # not_cell_lag <- identify_non_cells(el_lag, lagos)
        # el_lag[not_cell_lag]<-NA
        # 
        # 
        # 
        # el_del1 <- raster("1_input/data/elevation/delhi/ASTGTMV003_N28E076_dem.tif")
        # el_del2 <- raster("1_input/data/elevation/delhi/ASTGTMV003_N28E077_dem.tif")
        # el_del <- raster::merge(el_del1, el_del2)
        # el_del <- raster::crop(el_del, extent(delhi))
        # rm(el_del1, el_del2)
        # not_cell_del <- identify_non_cells(el_del, delhi)
        # el_del[not_cell_del]<-NA
        # 
        # 
        # 
        # 
        # 
         # 
        # save(w_rio, w_lag, w_del, pm_rio,pm_lag, pm_del,el_rio, el_lag, el_kat,el_del, file = "2_analysis/data/fig4-b-plotdata.RData" )
        # 
  
  
########################### plot ######################
        
        
  #load plot data
  load("data/fig4-b-plotdata.RData")
  
  #define color pals  
  pal_pm <- colorRampPalette(colors = MetBrewer::met.brewer("VanGogh2", 8))(2000) %>% rev()
  pal_elev <- colorRampPalette(colors = MetBrewer::met.brewer("Tiepolo", 8))(2000) %>% rev()
  pal_wealth <- colorRampPalette(colors = c("red3",MetBrewer::met.brewer("Benedictus", 16)[c(1:6, 8:13)], "navy"))(2000) 
  

  pdf("figures/fig4-b-maps.pdf", width = 12, height = 12)
        par(mfrow = c(3,3))
        par(mar = c(4,0,2,0))


        pm_rio_plot <- pm_rio
        pm_lag_plot <-  pm_lag

        pm_del_plot <-  pm_del

  ### RIO ###

        #[1]
        plot(rio, col = NA, border = NA)
        plot(pm_rio_plot, col = pal_pm, axes = F,add=T, legend = F)
        plot(rio, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "PM2.5",cex=2, line = -2)



        plotrix::gradient.rect(-43.72,-23.12,-43.21,-23.1,
                               col=(plotrix::smoothColors(colorRampPalette(pal_pm)(2000))),
                               gradient="x", border = 'black')




        #[2]

        int <- classIntervals(w_rio$rwi, style = "fixed", fixedBreaks = c(-10,seq(0.03, 1.7, .01),10))
        col_rio_w <- findColours(int, (pal_wealth))


        plot(rio, col = NA, border = NA)
        points(w_rio$longitude, w_rio$latitude, pch = 15, cex = 2.8, col = col_rio_w)
        plot(rio, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "Wealth",cex=2, line = -2)



        plotrix::gradient.rect(-43.72,-23.12,-43.21,-23.1,
                               col=(plotrix::smoothColors(colorRampPalette(pal_wealth)(2000))),
                               gradient="x", border = 'black')


        #[3]


        plot(rio, col = NA, border = NA)
        plot(el_rio, col = (pal_elev), axes = F,bty = "n",add=T, legend = F)
        plot(rio, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "Elevation",cex=2, line = -2)



        plotrix::gradient.rect(-43.72,-23.12,-43.21,-23.1,
                               col=(plotrix::smoothColors(colorRampPalette(pal_elev)(2000))),
                               gradient="x", border = 'black')








        ###### Delhi ########
        del <- delhi
        #[1]
        plot(del, col = NA, border = NA)
        plot(pm_del_plot, col = pal_pm, axes = F,bty = "n",add=T, legend = F)
        plot(del, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "PM2.5",cex=2, line = -2)

        plotrix::gradient.rect(85.22,27.535,85.52,27.55,
                               col=(plotrix::smoothColors(colorRampPalette(pal_pm)(2000))),
                               gradient="x", border = 'black')


        #[2]
        int <- classIntervals(w_del$rwi, style = "fixed", fixedBreaks = c(-10,seq(0.2, 1.6, .01),10))
        col_del_w <- findColours(int, (pal_wealth))

        plot(del, col = NA, border = NA)
        #plot(pm_del_plot, col = pal_pm, axes = F,bty = "n",add=T, legend = F)
        points(w_del$longitude, w_del$latitude, pch = 15, cex = 5, col = col_del_w)

        plot(del, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "Wealth",cex=2, line = -2)

        plotrix::gradient.rect(85.22,27.535,85.52,27.55,
                               col=(plotrix::smoothColors(colorRampPalette(pal_wealth)(2000))),
                               gradient="x", border = 'black')


        #[3]

        plot(del, col = NA, border = NA)
        plot(el_del, col = (pal_elev), axes = F,bty = "n",add=T, legend = F)
        plot(del, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "Elevation",cex=2, line = -2)

        plotrix::gradient.rect(85.22,27.535,85.52,27.55,
                               col=(plotrix::smoothColors(colorRampPalette(pal_elev)(2000))),
                               gradient="x", border = 'black')




####### LAGOS #########


      #[1]
        plot(lagos, col = NA, border = NA, ylim = c(6,7))
        plot(pm_lag_plot, col = pal_pm, axes = F,bty = "n",add=T, legend = F)
        plot(lagos, add = T, lwd =2)

          mtext(side = 3, adj = 0, text = "PM2.5",cex=2, line = -2)



        plotrix::gradient.rect(2.75,6.24,4.25,6.3,
                               col=(plotrix::smoothColors(colorRampPalette(pal_pm)(2000))),
                               gradient="x", border = 'black')



        #[2]

        int <- classIntervals(w_lag$rwi, style = "fixed", fixedBreaks = c(-10,seq(-.7, 1.6, .01),10))
        col_lag_w <- findColours(int, (pal_wealth))

        plot(lagos, col = NA, border = NA, ylim = c(6,7))
        #plot(pm_lag_plot, col = pal_pm, axes = F,bty = "n",add=T, legend = F)
        points(w_lag$longitude, w_lag$latitude, pch = 15, cex = 2, col = col_lag_w)

        plot(lagos, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "Wealth",cex=2, line = -2)



        plotrix::gradient.rect(2.75,6.24,4.25,6.3,
                               col=(plotrix::smoothColors(colorRampPalette(pal_wealth)(2000))),
                               gradient="x", border = 'black')


        #[3]




        plot(lagos, col = NA, border = NA, ylim = c(6,7))
        plot(el_lag, col = (pal_elev), axes = F,bty = "n",add=T, legend = F)
        plot(lagos, add = T, lwd =2)
        mtext(side = 3, adj = 0, text = "Elevation",cex=2, line = -2)



        plotrix::gradient.rect(2.75,6.24,4.25,6.3,
                               col=(plotrix::smoothColors(colorRampPalette(pal_elev)(2000))),
                               gradient="x", border = 'black')



        dev.off()


        
  
        
  ########## scatters ###########
        
        
        
        #PM vs wealth
        
        #rio
        rloc_rio <- w_rio[,c("quadkey","latitude","longitude","rwi")]
        pmloc_rio <-  data.frame(coordinates(pm_rio), pm_rio[])
        names(pmloc_rio)<-c("longitude","latitude","pm")
        pmloc_rio <- pmloc_rio %>% drop_na(pm)

        nn_rio <- FNN::knn( train = pmloc_rio[,c("longitude","latitude")],
                  test = rloc_rio[,c("longitude","latitude")],
                  cl = 1:nrow(pmloc_rio),k = 1
                    )

        rloc_rio[,c("lon_pm","lat_pm")] <- pmloc_rio[attr(nn_rio, "nn.index"),c("longitude","latitude")]
        rloc_rio$cell_pm <- cellFromXY(pm_rio,as.matrix(rloc_rio[,c("lon_pm","lat_pm")] ))
        rloc_rio$pm <- pm_rio[rloc_rio$cell_pm]
        
        

        scatter_pm_wealth_rio <- data.frame(wealth = rloc_rio$rwi, pm = rloc_rio$pm) %>% drop_na(pm, wealth)

        
        #delhi
        rloc_del <- w_del[,c("quadkey","latitude","longitude","rwi")]
        pmloc_del <-  data.frame(coordinates(pm_del), pm_del[])
        names(pmloc_del)<-c("longitude","latitude","pm")
        pmloc_del <- pmloc_del %>% drop_na(pm)
        
        nn_del <- FNN::knn( train = pmloc_del[,c("longitude","latitude")],
                            test = rloc_del[,c("longitude","latitude")],
                            cl = 1:nrow(pmloc_del),k = 1
        )
        
        rloc_del[,c("lon_pm","lat_pm")] <- pmloc_del[attr(nn_del, "nn.index"),c("longitude","latitude")]
        rloc_del$cell_pm <- cellFromXY(pm_del,as.matrix(rloc_del[,c("lon_pm","lat_pm")] ))
        rloc_del$pm <- pm_del[rloc_del$cell_pm]
        
        
        scatter_pm_wealth_del <- data.frame(wealth = rloc_del$rwi, pm =  rloc_del$pm) %>% drop_na(pm, wealth)
        
        
        
        #Lagos
        rloc_lag <- w_lag[,c("quadkey","latitude","longitude","rwi")]
        pmloc_lag <-  data.frame(coordinates(pm_lag), pm_lag[])
        names(pmloc_lag)<-c("longitude","latitude","pm")
        pmloc_lag <- pmloc_lag %>% drop_na(pm)
        
        nn_lag <- FNN::knn( train = pmloc_lag[,c("longitude","latitude")],
                            test = rloc_lag[,c("longitude","latitude")],
                            cl = 1:nrow(pmloc_lag),k = 1
        )
        
        rloc_lag[,c("lon_pm","lat_pm")] <- pmloc_lag[attr(nn_lag, "nn.index"),c("longitude","latitude")]
        rloc_lag$cell_pm <- cellFromXY(pm_lag,as.matrix(rloc_lag[,c("lon_pm","lat_pm")] ))
        rloc_lag$pm <- pm_lag[rloc_lag$cell_pm]
        
        
        scatter_pm_wealth_lag <- data.frame(wealth = rloc_lag$rwi, pm =  rloc_lag$pm) %>% drop_na(pm, wealth)
        
        
    ### elevation vs pm
        
        el_rio_pm <- resample(el_rio, pm_rio)
        scatter_pm_el_rio <- data.frame( cell_pm = 1:ncell(el_rio_pm), elev = el_rio_pm[], pm = pm_rio[]) %>% drop_na(pm, elev)


        el_lag_pm <- resample(el_lag, pm_lag)
        scatter_pm_el_lag <- data.frame( cell_pm = 1:ncell(el_lag_pm), elev = el_lag_pm[], pm = pm_lag[]) %>% drop_na(pm, elev)
        
        el_del_pm <- resample(el_del, pm_del)
        scatter_pm_el_del <- data.frame( cell_pm = 1:ncell(el_del_pm), elev = el_del_pm[], pm = pm_del[]) %>% drop_na(pm, elev)
        
        
        
        rloc_rio <- left_join(rloc_rio, scatter_pm_el_rio[,c("cell_pm","elev")])
        rloc_del <- left_join(rloc_del, scatter_pm_el_del[,c("cell_pm","elev")])
        rloc_lag <- left_join(rloc_lag, scatter_pm_el_lag[,c("cell_pm","elev")])
        
        
        
        
        
        
        
        #fig 4 (note also generates pm vs elevation scatters that didn't end up being used in final version of fig)
        pdf("figures/fig4-b-scatters.pdf", width = 12, height = 12)                
        
        par(mfrow =c(3,3))
        
        plot(scatter_pm_wealth_rio, ylim = c(10,20), xlim = c(-0.5, 2), axes = F, xlab = "", ylab = "", col = NA)
        points(scatter_pm_wealth_rio, pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x  =seq(-0.5, 2, .01), y =  predict(loess(pm ~ wealth, data = scatter_pm_wealth_rio), newdata = data.frame(wealth = seq(-0.5, 2, .01))), col  = 'blue',lwd=3)
        mtext(side = 2, text = "PM2.5",line=3,cex=2)
        mtext(side = 1, text = "Wealth",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
    
        
        
        
        plot(scatter_pm_wealth_del, axes = F, xlab = "", ylab = "", col = NA)
        points(scatter_pm_wealth_del, pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x  =seq(-0.4, 1.8, .01), y =  predict(loess(pm ~ wealth, data = scatter_pm_wealth_del), newdata = data.frame(wealth = seq(-0.4, 1.8, .01))), col  = 'blue',lwd=3)
        mtext(side = 2, text = "PM2.5",line=3,cex=2)
        mtext(side = 1, text = "Wealth",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
        plot(scatter_pm_wealth_lag, axes = F, xlab = "", ylab = "", col = NA)
        points(scatter_pm_wealth_lag, pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x  =seq(-1, 1.9, .01), y =  predict(loess(pm ~ wealth, data = scatter_pm_wealth_lag), newdata = data.frame(wealth = seq(-1, 1.9, .01))), col  = 'blue',lwd=3)
        mtext(side = 2, text = "PM2.5",line=3,cex=2)
        mtext(side = 1, text = "Wealth",line=3,cex=2)
        #text(x = -0.75, y = 39, labels = paste("slope = ",round(coef(lm(pm ~ wealth, data = scatter_pm_wealth_lag))[2],1), sep =""),cex=1.5)
        axis(1)
        axis(2, las=2)
        
        
        ###
        
        plot(scatter_pm_el_rio[,2:3], ylim = c(10,20), xlim = c(0,700), axes = F, xlab = "", ylab = "", col = NA)
        points(scatter_pm_el_rio[,2:3], pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x = 0:700, y = predict(loess(pm ~ elev, data = scatter_pm_el_rio), newdata= data.frame(elev = 0:700)), col  = 'blue',lwd=3)
        mtext(side = 2, text = "PM2.5",line=3,cex=2)
        mtext(side = 1, text = "Elevation",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
        
        plot(scatter_pm_el_del[,2:3], axes = F, xlab = "", ylab = "", col = NA)
        points(scatter_pm_el_del[,2:3], pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x = 175:300, y = predict(loess(pm ~ elev, data = scatter_pm_el_del), newdata= data.frame(elev = 175:300)), col  = 'blue',lwd=3)
        mtext(side = 2, text = "PM2.5",line=3,cex=2)
        mtext(side = 1, text = "Elevation",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
        
        plot(scatter_pm_el_lag[,2:3], axes = F, xlab = "", ylab = "", col = NA)
        points(scatter_pm_el_lag[,2:3], pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x = 0:700, y = predict(loess(pm ~ elev, data = scatter_pm_el_lag), newdata= data.frame(elev = 0:700)), col  = 'blue',lwd=3)
        mtext(side = 2, text = "PM2.5",line=3,cex=2)
        mtext(side = 1, text = "Elevation",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
        
        
        
        
        
        plot(rloc_rio[,c("rwi","elev")], ylim = c(0,600), xlim = c(-0.6,1.9), axes = F, xlab = "", ylab = "", col = NA)
        points(rloc_rio[,c("rwi","elev")], pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x = seq(-0.5, 2, .01), y = predict(loess(elev ~ rwi, data = rloc_rio), newdata= data.frame(rwi = seq(-0.5, 2, .01))), col  = 'blue',lwd=3)
        mtext(side = 1, text = "Elevation",line=3,cex=2)
        mtext(side = 2, text = "Wealth",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
        
        plot(rloc_del[,c("rwi","elev")], axes = F, xlab = "", ylab = "", col = NA)
        points(rloc_del[,c("rwi","elev")], pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x = seq(-0.5, 2, .01), y = predict(loess(elev ~ rwi, data = rloc_del), newdata= data.frame(rwi = seq(-0.5, 2, .01))), col  = 'blue',lwd=3)
        mtext(side = 1, text = "Wealth",line=3,cex=2)
        mtext(side = 2, text = "Elevation",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
        
        plot(rloc_lag[,c("rwi","elev")], axes = F, xlab = "", ylab = "", col = NA)
        points(rloc_lag[,c("rwi","elev")], pch = 21, bg = add.alpha('gray70', .9), col = add.alpha('white', 1), cex=2, lwd = 0.25)
        lines(x = seq(-0.5, 2, .01), y = predict(loess(elev ~ rwi, data = rloc_lag), newdata= data.frame(rwi = seq(-0.5, 2, .01))), col  = 'blue',lwd=3)
        mtext(side = 1, text = "Wealth",line=3,cex=2)
        mtext(side = 2, text = "Elevation",line=3,cex=2)
        axis(1)
        axis(2, las=2)
        
        
        
        dev.off()
        
        
        
        
        
        
        
